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1.  INTRODUCTION 


» 


This  report  presents  an  algorithm  for  computing  the  stationary- 
point  of  a quadratic  function  of  n variables  subject  to  a set  of  m(m  < n) 
linear  equality  constraints.  The  procedure  is  compact  in  the  sense  that  it 
requires  no  two-dimensional  arrays  of  computer  storage  beyond  that  needed 
to  store  the  problem  data.  The  use  of  a Householder  orthogonal  decom- 
position by  this  method  should  not  degrade  the  numerical  conditioning  of  the 
original  problem.  This  method  is  applicable  to  problems  with  singular  < 

Hessian  matrices,  and  can  be  adapted  for  use  in  a general  quadratic 
programming  algorithm. 

In  the  subsequent  sections  of  this  report,  the  identifying  numbers 
of  equations  in  the  text  are  enclosed  with  parentheses,  and  the  identifying 
numbers  of  references  are  enclosed  with  brackets. 

I 
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When  A is  positive  (negative)  definite,  the  solution  defines  the  unique 
stationary  point  which  corresponds  to  the  minimum  (maximum)  of  f restricted 
to  the  constraint  surface.  While  the  stated  problem  may  be  of  interest  by 
itself,  typically  it  may  appear  as  a subproblem  in  a more  general  application. 
For  example,  many  quadratic  programming  algorithms  solve  a series  of 
problems  of  this  type  with  different  constraint  sets.  Furthermore,  an 
algorithm  designed  to  optimize  a non-quadratic  function  subject  to  nonlinear 
constraints  may  pose  a series  of  quadratic -linear  problems  to  approximate 
the  behavior  of  the  actual  functions.  Consequently,  it  is  desirable  to  develop 
a computational  algorithm  which  will  compute  a solution  to  the  problem  with- 
out restricting  the  rank  of  A.  The  computed  solution  should  be  the  unique 
solution  when  it  exists,  and  should  be  uniquely  defined  by  the  algorithm  when 
a unique  solution  does  not  exist. 

When  the  minimum  norm  of  the  projected  gradient  is  zero,  the  solution 
to  the  stated  problem  is  a stationary  point  of  the  Lagrangian  function 


t 


(Cx  - d) 


(3) 


L(x,  X)  = Ax  + b^*  x + XT 

where  X.  is  the  m-vector  of  Lagrange  multipliers.  Setting  the  derivative  of 
this  function  with  respect  to  x and  X equal  to  zero  yields  the  set  of  linear 
equations 


The  optimal  solution  x and  the  corresponding  multiplier  values  X 

can  be  obtained  by  solving  the  system  (4).  It  is  not  necessary  to  assume  that 

A is  of  full  rank.  When  the  problem  has  a unique  solution,  the  system  may 

be  solved  using  a suitable  algorithm  for  linear  equations.  If  the  possibility 

of  a non-unique  solution  exists,  the  system  may  be  solved  as  a linear  least 

1 2 

squares  problem.  This  approach  has  been  used  ’ utilizing  the  linear 

3 

least  squares  algorithm  of  Hanson  and  Lawson  . A defect  in  this  approach 
is  the  need  to  store  the  (n  + m)  x (n  + m)  coefficient  array. 

An  approach  which  can  be  implemented  using  only  the  storage  required 
for  A and  C can  be  derived  by  inverting  the  coefficient  matrix  in  a partitioned 


*J.  T.  Betts,  "A  Gradient  Projection  - Multiplier  Method  for  Nonlinear 
Programming,  " Journal  of  Optimization  Theory  and  Applications, 

Vol  24,  No.  4 (April  1978). 

2 

J.  T.  Betts,  "An  Accelerated  Multiplier  Method  for  Nonlinear 
Programming,  " Journal  of  Optimization  Theory  and  Applications, 

Vol  21,  No.  2 (February  1977). 
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R.  J.  Hanson  and  C.  L.  Lawson,  "Solving  Least  Squares  Problems,  " 
Prentice -Hall,  Englewood  Cliffs,  N.  J.  (1974). 
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form.  It  is  easily  demonstrated  that 


'a 

cT- 

— A 

'B1 

V 

c 

o 

= 

b,T 

B, 

L 2 

3J 

where 

-1  -I  T -1  -1 

Bj  = A 1 - A C1  M 1 CA  1 

-1  T -1 
B2  = A C M 

B3  = M"1 

and 

-1  T 

M = CA  C 


(5) 


(6) 

(7) 

(8) 


(9) 


The  partitioned  form  of  the  inverse  plays  an  important  role  in  a number  of 
quadratic  programming  algorithms,  including  those  of  Goldfarb4  and 
Fletcher  , as  well  as  in  the  constrained  minimization  algorithm  of 
Murtagh  and  Sargent  . Two  significant  points  deserve  comment  regarding 
this  approach.  First,  if  it  is  assumed  that  A-1  exists,  the  submatrices  in  (5) 
can  be  computed  directly.  Goldfarb,  Murtagh,  and  Sargent  assume  that  A is 


4 

D.  Goldfarb,  "Extension  of  Newton's  Method  and  Simplex  Methods  for 
Solving  Quadratic  Programs,  " In  Numerical  Methods  for  Nonlinear 
Optimization,  F.  A.  Lootsman  (Ed.)  Academic  Press,  London,  ch.  17 
(1972). 

5 

R.  Fletcher,  "A  General  Quadratic  Programming  Algorithm,  " 

J.  Inst.  Math.  Appl. , Vol  8 (1971). 

^B.  A.  Murtagh  and  R.  W.  H.  Sargent,  "A  Constrained  Minimisation 
Method  for  Quadratic  Convergence,  " In  Optimisation,  R.  Fletcher 
(Ed.),  Academic  Press,  London,  ch.  14  (1969) 


positive  definite,  which  ensures  that  A is  of  full  rank.  If  the  algorithm  is  part 
of  a more  general  nonlinear  programming  algorithm  as  in  [1]  and  [2],  it  is 
restrictive  to  assume  that  A is  of  full  rank.  For  example,  the  approach  could 
not  be  used  as  part  of  an  algorithm  to  optimize  a linear  function  subject  to 
nonlinear  constraints.  In  these  general  applications,  it  is  really  only 
necessary  that  the  Hessian  matrix  restricted  to  the  constraint  surface  be 
definite.  Fletcher  does  not  assume  A is  definite,  noting  that  the  partitions 
Bl'  B2*  and  B3  must  exist  if  the  solution  to  the  problem  is  unique.  However, 
to  compute  the  initial  submatrices  in  the  computer  implementation  of  his 
quadratic  programming  algorithm  , it  is  necessary  to  invert  the  full 
(n  + m)  x (n  + m)  matrix. 

Even  if  A is  assumed  to  be  definite,  the  partitioned  approach  to  the 

T 

problem  suffers  from  a second  defect.  This  occurs  when  A = I,  M = CC  , 
which  is  referred  to  as  the  normal  matrix.  In  this  case,  the  condition 
number  of  M is  the  square  of  the  condition  number  of  C and  it  is  generally 
recognized  that  the  formation  of  M is  to  be  avoided. 

In  summary,  direct  solution  of  the  system  (4)  is  not  compact  from  a 
storage  standpoint.  The  various  forms  of  solving  the  partitioned  system  (5), 
although  compact,  require  operations  which  can  degrade  the  numerical 
conditioning  of  the  given  problem  and  are  arbitrarily  restrictive  with  regard 
to  Hessian  matrix  A.  A new  algorithm  will  be  proposed  which  is  compact, 
does  not  degrade  the  numerical  conditioning,  and  makes  no  restrictions 
concerning  the  rank  of  A. 


7 

R.  Fletcher,  "A  FORTRAN  Subroutine  for  Quadratic  Programming, " 
Report  No.  AERE  R6370,  UKAEA  Research  Group,  Harwell,  England. 
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3.  ORTHOGONAL  DECOMPOSITION  ALGORITHM 

In  this  section,  an  algorithm  is  developed  for  solving  the  stated 
constrained  stationary  point  problem  using  an  orthogonal  decomposition  of 
the  constraint  matrix.  The  algorithm  is  an  extension  of  the  linearly  con- 
strained linear  least  squares  algorithm  LSE  given  in  [3],  and  makes  use 
of  Theorem  (3,  19)  and  Theorem  (2.3)  stated  therein. 

Define  the  orthogonal  decomposition  of  C by 

C = RKT  (10) 

where  K is  an  n x n orthogonal  matrix  and  R 

R = [Rn  O]  (11) 

where  Rjjisanmxm  nonsingular  triangular  matrix.  Substituting  (10) 
into  (2) 

RKT  x = d.  (12) 

Define  the  n-vector  y by 

y = KTx  (13) 


and  the  partitions  of  K and  y 


K = [Kj  K2] 
m n-m 


(14) 


(15) 


From  (11),  (13),  and  (15),  one  can  write  (12)  as 


RK  x = Ry  = [Rn  O]  y*j  = R^j  = d (16) 

Since  Rjj  i **  non-singular,  (16)  determines  the  m-vector  y^.  Call  the  solution 
y^.  The  (n-m)  -vector  y2  is  arbitrary. 

Pre -multiply  (13)  by  K to  give 

Ky  = KKTx  = x (17) 

T 

where  KK  = I,  since  K is  orthogonal.  From  (14)  and  (15),  (17)  becomes 

X = Ky  = [Kjiy  y*  = KjYi  + K2y2  (18) 

Observe  that  all  points  satisfying  the  constraints  in  Eq.  (2)  can  be  represented  as 
functions  of  the  n-m  parameters  y2  when  the  solution  of  (16)  y^  is  substituted 
in  Eq.  (18).  In  fact  if  there  were  no  other  conditions  to  satisfy,  a reasonable 
choice  for  the  arbitrary  parameters  y2  would  be  zero,  in  which  case  x is  the 
minimum  norm  solution  to  the  constraints. 

Instead  of  setting  y2  = 0,  the  choice  of  y2  shall  be  determined  by  a 
different  criterion.  Substitute  the  parametric  representation  of  x from  (18) 
with  y^  = yj  into  (1)  to  obtain 


*■  = * (KjVj  + K2y2)T  A (KLy  + K^) 


+ b (Vi + K2y2)  • 


The  gradient  with  respect  to  the  variables  y2  is 

VI  = K2TA(K1^1  + K2y2)  + K2Tb. 
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Let  us  define  to  be  the  value  of  which  minimizes 

11**11  = ||K2TAK2y2  + (K2Tb  + K^AK^U  (21) 

and  is  of  minimum  norm,  i.  e. , minimizes  ||  y2  If  a stationary  point  of  f 
restricted  to  the  constraint  surface  exists,  then  ||  7f  ||  = 0 and  y defines  the 
optimal  value.  If  the  matrix  K 2 AK2  is  indefinite,  the  minimum  norm 
criterion  uniquely  determines  y_.  In  fact,  when  the  rank  of  K^T  AK,  is  zero 
as  is  the  case  for  a linear  objective  function,  the  solution  which  minimizes 
the  norm  of  ||  y2  ||  is  just  y2  = 0.  Notice  also  that  the  formation  of  the 
matrix  K2T  A K2  should  not  degrade  the  conditioning  of  the  original  problem. 
Observe  also  that  a solution  which  minimizes  ||  y2||  minimizes  II  ^2y2  11  since 

II  Y2  II  = II  k2V2  N*  and  K2y2  iS  ^USt  the  orthog°nal  component  of  x in  the 
constraint  surface. 

In  summary,  the  original  constrained  optimization  problem  is  replaced 
by  a lower  dimensional  unconstrained  least  squares  problem  in  the  variables 
y2,  after  choosing  the  variables  yj  to  satisfy  the  constraints.  The  method 
has  the  property  that  the  unique  minimum  length  solution  of  the  derived 
unconstrained  problem  defines  the  unique  solution  of  the  constrained  problem 
when  it  exists.  When  the  constrained  problem  has  no  unique  solution,  the 
algorithm  computes  a unique  point  which  satisfies  the  constraints,  minimizes 
the  norm  of  the  gradient  on  the  constraint  surface,  and  minimizes  the  length 
of  the  orthogonal  component  in  the  constraint  surface. 


4.  COMPUTATIONAL  ALGORITHM 


In  this  section  a computational  procedure  based  on  the  approach  derived 
in  Section  3 is  developed.  The  procedure  is  organized  so  that  no  additional 
two-dimensional  arrays  are  needed.  Specifically,  the  original  problem  data 
stored  in  A,  C and  b is  destroyed  by  the  algorithm.  Quantities  written  with 
a tilde  can  replace  quantities  without  a tilde  in  storage,  and  quantities  written 
with  a circumflex  can  overwrite  quantities  written  with  a tilde. 

Step  1.  Compute  the  orthogonal  matrix  K and  postmultiply  C by  it  to 
triangularize  C,  i.  e. , 


CK  : 

r ** 

= r°i 

O] 

m 

(22) 

m n-m 

Step  2. 

Compute 

^0 

A = 

T 

K A 

(23) 

Step  3. 

Form  the  last  n-m  rows 

of  the  matrix 

A 

A where 

A 

A = 

0*0 

AK 

(24) 

Observe  that  from  (23)  and  (24) 

Vak, 

KlTAK2' 

1 m 

A 

A = 

KT  AK 

= 

[k2T  AKj 

K2^K2- 

J n-m 

~rn 

n-m 

(25) 

A11 

A12 

A 

La2i 

A 

A22 

Step  4.  Compute 


^ T 

b = K d 


Step  5.  Solve  the  lower  triangular  system 

eiyl  = d 

for  the  m-vector  y^. 

Step  6.  Compute 


-b2  - A21yl 


where 


b = 


J2J! 


m 


i n-m 


S = 


Lbn 


m 


n-m 


Step  7.  Determine  y^  as  the  minimum  length  solution  of  the 
least  square  problem 


min 


22  y2 


Observe  that  this  process  is  equivalent  to  solving  (21) 
Step  8.  Construct  the  solution  vector 


x = K y 


using  yj  as  computed  in  Step  5,  and  y^  from  Step  7. 
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(26) 


(27) 


(28) 


(29) 

linear 


(30) 


(31) 


The  algorithm  described  has  been  implemented  in  the  subroutine 
HSQP.  The  FORTRAN  listing  of  this  subroutine  is  found  in  [8],  The  sub- 
routine makes  extensive  use  of  the  subroutines  HFTI  and  H12  which 
implement  the  algorithms  referred  to  as  HFTI,  HI,  and  H2  in  [3],  The 
subroutine  HFTI  computes  the  minimum  length  solution  to  a linear  least 
squares  problem.  HFTI  requires  storage  for  the  problem  data  and  three 
one -dimensional  work  arrays.  The  subroutine  H12  implements  algorithm 
HI  and  H2  for  the  construction  and  application  of  a Householder  transform- 
ation. Using  H12  it  is  not  necessary  to  explicitly  form  the  orthogonal  matrix 
K of  Eq.  (10).  Instead,  the  elements  necessary  to  construct  the  matrix  can 
be  stored  in  the  upper  triangular  portion  of  the  original  matrix  C and  some 
one-dimensional  work  arrays.  Successive  applications  of  the  matrix  K to 
other  vectors  implicitly  reconstruct  the  original  transformations. 

The  total  storage  required  for  subroutine  HSQP,  including  that  re- 
quired to  specify  the  problem  data,  is  = n(m  + n)  + 5n  - m + 4.  In  contrast, 
any  algorithm  which  solves  (4)  directly  will  require  at  least  = (n  + m) 

(n  + m + 1)  storage  locations.  Consequently,  for  some  problems  can  be 
nearly  twice  as  large  as  Nj,  The  algorithm  is  used  repeatedly  as  part  of 
the  general  nonlinear  programming  algorithm  described  in  [1],  In  particular, 
all  of  the  extrapolation  steps  used  in  the  constraint  phase  of  this  algorithm 
employ  HSQP.  Computational  experience  with  the  algorithm  includes  its 
use  to  solve  the  set  of  17  equality  constrained  and  34  inequality  constrained 
problems  in  [1],  as  well  as  a number  of  larger  engineering  applications. 

One  typical  application  is  described  in  [9],  The  largest  engineering  applica- 
tion of  the  algorithm  to  date  occurred  in  an  optimum  solid  rocket  motor 
design  problem  which  involved  48  variables  and  83  constraints. 


J.  T.  Betts,  "Algorithm  (To  Appear):  The  Stationary  Point  of  a Quadratic 
Function  Subject  to  Linear  Constraints,  " ACM  Trans.  Math.  Software. 

9 11 

J.  T.  Betts,  "Optimal  Three  Burn  Orbit  Transfer,  " AIAA  Journal. 

Vol  15,  No.  6 (June  1977),  


-15- 


5.  SUMMARY 

An  algorithm  for  computing  the  stationary  point  of  a quadratic  function 
of  n variables  subject  to  m linear  equality  constraints  is  developed.  The 
algorithm  has  been  implemented  in  FORTRAN.  The  implementation  is 
compact  since  it  requires  no  two-dimensional  arrays  beyond  that  needed  to 
define  the  problem.  The  algorithm  avoids  mathematical  operations  which 
would  degrade  the  conditioning  of  the  original  problem  by  utilizing  an 
orthogonal  decomposition  of  the  constraint  matrix.  The  solution  generated 
by  the  algorithm  is  characterized  by  three  properties:  (a)  the  constraints 
are  satisfied,  (b)  the  norm  of  the  gradient  of  the  objective  function  restricted 
to  the  constraint  surface  is  minimized  and,  (c)  among  all  solutions  satisfying 
the  first  two  properties,  the  minimum  length  solution  is  chosen.  When  the 
stated  problem  has  a unique  solution,  satisfaction  of  the  first  two  properties 
defines  the  point.  Nevertheless,  the  algorithm  is  not  restricted  to  problems 
with  definite  Hessian  matrices.  The  algorithm  has  been  successfully  tested 
as  part  of  a general  nonlinear  programming  algorithm. 
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APPENDIX 

THE  STATIONARY  POINT  OF  A QUADRATIC  FUNCTION 
SUBJECT  TO  LINEAR  CONSTRAINTS 


This  algorithm  implements  the  method  developed  in  the  preceding 
sections  of  this  report. 


SUBROUTINE  HS'JP  (A,B,C,D,M,  N,TAU,G,  H,U,  IP, H AX R A, MAXRC , D JNORM , X, 
$ KRANK) 


DIMENSION  B(1)  , D ( 1 ) , G ( 1 ) ,H(1)  ,0(1)  , IP  (1)  ,DJN0RM  (1)  ,X  (1) 

DIMENSION  A (HAXRA, 1) ,C (MAXRC, 1) 

C PROGRAMMER  AND  DATE:  J.T. BETTS,  JAN.  1978. 

C 

C PURPOSE:  GIVEN  AN  M X N MATRIX  C (OF  RANK  M) , AN  M VECTOR  D, 

C AN  N X N SYMMETRIC  MATRIX  A,  AND  AN  N VECTOR  B,  FIND  THE 

C STATIONARY  POINT  X OF  THE  (JOADRATIC 

C 

C J * .5*  (X**T)  *A*X  ♦ (B**T)*X 

C 

C SUBJECT  TO  THE  CONSTRAINTS 

C 

C C*X  = D. 

C 

C IF  A STATIONARY  POINT  DOES  NCT  EXIST  THE  ALGORITHM  WILL  FIND 

C A POINT  WHICH  SATISFIES  THE  CONSTRAINTS  AND  MINIMIZES  THE 

C NORM  OF  THE  GRADIENT  OF  J PROJECTED  ON  THE  CONSTRAINT  SURFACE. 

C 

C ALGORITHM:  ORTHOGONAL  DECOMPOSITION  OF  C MATRIX  USING 

C HOUSEHOLDER  TRANSFORMATIONS,  FOLLOWED  BY  APPLICATION  OF  THE 


non 
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C OPTIMALITY  CONDITIONS  IN  THE  FfDUCED  VARIABLES. 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


INPUT: 

A N X N SYMMETRIC  HESSIAN  MATRIX 

B N DIMENSIONAL  GRADIENT  VECTOR 

C M X N JACOBIAN  MATRIX  (RANK  M) 

D M DIMENSIONAL  CONSTRAINT  VECTOR 

M THE  NUMBER  OF  CONSTRAINTS 

N THE  NUMBER  OF  VARIABLES 

TAU  PSEUDORANK  TEST  PARAMETER.  FOB  A MACHINE  WITH  K 

SIGNIFICANT  FIGURES  AN  APPROPRIATE  VALUE  IS 
TAU  = 1. E- (K+2) . 

G AUXILLIARY  STORAGE  (LENGTH  M) 

H AUXILLIARY  STORAGE  (LENGTH  N-H) 

U AUXILLIARY  STORAGE  (LENGTH  N-H) 

IP  AUXILLIARY  STORAGE  (LENGTH  N-M) 

MAX  RA  MAXIMUM  ROW  DIMENSION  OF  A ( MAXRA  N) 

MAXRC  MAXIMUM  ROW  DIMENSION  OF  C (HAXRC  M) 

OUTPUT: 

DJNORM  PROJECTED  GRADIENT  NORM  (EERO  IF  X IS  A STATIONARY 
POINT,  NEGATIVE  IF  THERE  IS  AN  INPUT  ERROR) 

X COMPUTED  STATIONARY  POINT 

KRANK  PS5UDOPANK  OF  PROJECTED  HESSIAN  MATRIX  (K2**T) *A*K2. 

WHEN  KRANK  .LT.  N-M  THE  PROJECTION  OF  X ON  THE 
CONSTRAINT  SURFACE  HAS  MINIMUM  NORM. 

NOTE:  THE  INPUT  VALUES  OF  A,B,C,  AND  D ARE  DESTROYED. 


INITIALIZATION 

KRANK  = 0 
MP1  = M ♦ 1 
NMM  = N - M 
DJNORM  (1 ) = -1. 


CHECK  FOR  INPUT  ERRORS 


C 


IF  (N.  EQ.  O.OR.N.GT.  MA  XRA  . OR  . M.  GT.  MAXRC.  0 R.  M . GT.  N)  RETURN 
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IF  THE 
IF  (M.  EQ.  0) 


PROBLEM  IS  UNCONSTRAINED  GO  TO  STEP  7 
GO  TO  100 


STEP.  1.  COMPUTE  ORTHOGONAL  MATRIX  K.  TRIANGULARIZE  C. 


10 


DO  10  I = 1 r M 

CALL  H12  (1  ,I,I*1,N,C  (I,  1)  , MAXRC,G(I)  ,C  (1*1 ,1)  , MAXRC, 1 , M-I) 
CONTINUE 

IF  (H.  EQ.  N)  GO  TO  50 


STEP  2.  COMPOTE  ATILDA  = (K**T) *A 
DO  20  I = 1,N 

CALL  H12  (2,I,I*1,N,C  (I,  1)  ,HAXRC,G(I)  , A,  1.MAXRA.N) 

CONTINUE 


STEP  3.  FORM  THE  LAST  N-M  RONS  OF  AHAT  = ATILDA*K;  I.E. 
COMPOTE  A21HAT  * (K2**T) *A*K1  AND  A22HAT  * (K2**T) *A*K 2 

DO  30  I s 1,(1 

CALL  H12  (2,I,I*1,N,C(I,1)  .MAXRC, G(I)  ,A(MP1,1)  , MAXR  A , 1 , NMM) 

33  CONTINUE 


STEP  4.  COMPUTE  BTILDA  = <K»*T) *B 
DO  40  I = 1 , M 

CALL  H 12  (2,I,I»1,N,C(I,  1)  , HAXRC,G<I)  ,B,  1,1,1) 
40  CONTINUE 


STEP  5.  COMPOTE  XI  HAT  BT  SOIVING  THE  LONER  TRIANGULAR 
SYSTEM  C*X1  * D.  STORE  IN  X. 

50  CONTINUE 
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X(1)  = D (1 ) /C  ( 1 , 1) 

IF(M.EQ.  1)  GO  TO  90 
DO  70  I = 2,N 
III  1 * I - 1 
x CD  = D (I) 

DO  60  J = 1 , IN  1 

X(I)  = X (X)  - C(I,J)*X(J) 

60  CONTINUE 

x ( I)  = X(I)/C(I,I) 

70  CONTINUE 
90  CONTINUE 

WHEN  THERE  ARE  NO  DEGREES  OF  FREEDOM  GO  TO  STEP  8 
IF  (N.  EQ.  N)  GO  TO  140 


STEP  6.  COMPOTE  B2HAT  = -B2TTLDA  - A21HAT*Y1HAT 

DO  90  I = MP 1 , N 
3(1)  = -B  (I) 

DO  93  J * 1,11 

3(1)  = B (I ) - A (I»  0)  *X  (3) 

CONTINOE 


STEP  7.  SOLVE  A22HAT*Y2  = B2HAT  FOR  Y2  USING  HPTI 
100  CONTINOE 

COMPUTE  PSEUDORANK  TEST  PARAMETER  EPS 

EPS  = TAU 
DO  120  J = MP1,N 
COLNRM  * C. 

DO  110  I = MP1,N 
COLNRN  = COLNRN  ♦ A(I,J)**2 
110  CONTINUE 

EPS  * ANAX1 (EPS, TA0*SQBT (COLNRN) ) 

120  CONTINUE 

CALL  HFTI(A(NP1,NP1) ,RAXRA,NHH,NRH,B(HP1) , 1 , 1 , EPS, KR AN K, DJNORN , 


o I 


oonno  onn 
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$ H , U,  IP) 

c 

DO  130  I = NP1,N 
X(I)  = B (I) 

130  CONTINUE 

IF  THE  PROBLEM  IS  UNCONSTRAINED,  RETURN. 
IP(H.EQ.  C)  RETURN 


STEP  8.  COMPUTE  X = K*Y 

140  CONTINUE 

DO  150  K = 1,H 
I = M PI  - K 

CALL  H12  (2,I,I+1,N,C  (I,  1)  ,MAXRC,G(I)  ,X,  1,1,1) 
150  CONTINUE 
C 

RETURN 

END 
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